NASA/TM-20 14-2 1 8660 



Training Knowledge Bots for Physics-Based 
Simulations Using Artificial Neurals Networks 


Jay Ming Wong 

University of Massachusetts, Amherst, Massachusetts 
Jamshid A. Samareh 

Langley Research Center, Hampton, Virginia 


November 2014 




NASA STI Program ... in Profile 


Since its founding, NASA has been dedicated to the 
advancement of aeronautics and space science. The 
NASA scientific and technical information (STI) 
program plays a key part in helping NASA maintain 
this important role. 

The NASA STI program operates under the auspices 
of the Agency Chief Information Officer. It collects, 
organizes, provides for archiving, and disseminates 
NASA’s STI. The NASA STI program provides access 
to the NTRS Registered and its public interface, the 
NASA Technical Reports Server, thus providing one 
of the largest collections of aeronautical and space 
science STI in the world. Results are published in both 
non-NASA channels and by NASA in the NASA STI 
Report Series, which includes the following report 
types: 

• TECHNICAL PUBLICATION. Reports of 
completed research or a major significant phase of 
research that present the results of NASA 
Programs and include extensive data or theoretical 
analysis. Includes compilations of significant 
scientific and technical data and information 
deemed to be of continuing reference value. 

NASA counter-part of peer-reviewed formal 
professional papers but has less stringent 
limitations on manuscript length and extent of 
graphic presentations. 

• TECHNICAL MEMORANDUM. 

Scientific and technical findings that are 
preliminary or of specialized interest, 
e.g., quick release reports, working 

papers, and bibliographies that contain minimal 
annotation. Does not contain extensive analysis. 

• CONTRACTOR REPORT. Scientific and 
technical findings by NASA-sponsored 
contractors and grantees. 


• CONFERENCE PUBLICATION. 

Collected papers from scientific and technical 
conferences, symposia, seminars, or other 
meetings sponsored or 

co-sponsored by NASA. 

• SPECIAL PUBLICATION. Scientific, 
technical, or historical information from NASA 
programs, projects, and missions, often 
concerned with subjects having substantial 
public interest. 

• TECHNICAL TRANSLATION. 
English-language translations of foreign 
scientific and technical material pertinent to 
NASA’s mission. 

Specialized services also include organizing 
and publishing research results, distributing 
specialized research announcements and feeds, 
providing information desk and personal search 
support, and enabling data exchange services. 

For more information about the NASA STI program, 
see the following: 

• Access the NASA STI program home page at 
http://www.sti.nasa.gov 

• E-mail your question to help@sti.nasa. gov 

• Phone the NASA STI Information Desk at 
757-864-9658 

• Write to: 

NASA STI Information Desk 
Mail Stop 148 

NASA Langley Research Center 
Hampton, VA 23681-2199 



NASA/TM-20 14-2 1 8660 



Training Knowledge Bots for Physics-Based 
Simulations Using Artificial Neurals Networks 


Jay Ming Wong 

University of Massachusetts, Amherst, Massachusetts 
Jamshid A. Samareh 

Langley Research Center, Hampton, Virginia 


National Aeronautics and 
Space Administration 

Langley Research Center 
Hampton, Virginia 23681-2199 


November 2014 




The use of trademarks or names of manufacturers in this report is for accurate reporting and does not constitute an 
official endorsement, either expressed or implied, of such products or manufacturers by the National Aeronautics and 
Space Administration. 


Available from: 


NASA STI Program / Mail Stop 148 
NASA Langley Research Center 
Hampton, VA 23681-2199 
Fax: 757-864-6500 




1 


Training Knowledge Bots for Physics-Based 
Simulations Using Artificial Neural Networks 

Jay Ming Wong and Jamshid A. Samareh 
NASA Langley Research Center 
Hampton, Virginia 23681-2199 


Abstract — Millions of complex physics-based simulations are 
required for design of an aerospace vehicle. These simulations 
are usually performed by highly trained and skilled analysts, who 
execute, monitor, and steer each simulation. Analysts rely heavily 
on their broad experience that may have taken 20-30 years to 
accumulate. In addition, the simulation software is complex in na- 
ture, requiring significant computational resources. Simulations 
of system of systems become even more complex and are beyond 
human capacity to effectively learn their behavior. IBM has 
developed machines that can learn and compete successfully with 
a chess grandmaster and most successful jeopardy contestants. 
These machines are capable of learning some complex problems 
much faster than humans can learn. 

In this paper, we propose using artificial neural network to 
train knowledge bots to identify the idiosyncrasies of simulation 
software and recognize patterns that can lead to successful 
simulations. We examine the use of knowledge bots for applica- 
tions of computational fluid dynamics (CFD), trajectory analysis, 
commercial finite-element analysis software, and slosh propellant 
dynamics. We will show that machine learning algorithms can 
be used to learn the idiosyncrasies of computational simulations 
and identify regions of instability without including any additional 
information about their mathematical form or applied discretiza- 
tion approaches. 


I. Introduction 

Traditionally, individuals perform millions of simulations 
during the course of a design process by providing inputs, 
analyzing the simulation outputs, and identifying patterns 
that lead to successful simulations. Trained and experience 
analysts have learned the underlying simulation idiosyncrasies, 
and they are able to predict the simulation behavior given 
a set of input parameters, prior to actually performing the 
simulation. Humans are great at identifying complex patterns 
for events with small number of variables and given sufficient 
time to process the data and develop an expertise. Machine 
learning algorithms, such as artificial neural networks, are 
rapidly catching up and exceeding human capacity to iden- 
tify these patterns that leads to knowledge creation. This is 
especially useful for complex simulation as the underlying 
implementation may be too difficult to decipher. However, to 
obtain such knowledge requires years of simulation experience. 
Furthermore, the software itself may not only be extremely 
complex, but may require a long period of time for execution. 
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We propose an approach for a trained classifier, in which we 
refer to as a knowledge bot that learns the underlying properties 
of the black-box, much like a trained expert, and provide a 
quick response to whether or not a set of input parameters 
will result in failure before actually executing the black-box 
with inputs. Therefore, designers may be able to produce initial 
designs given the information of the knowledge bot , meanwhile 
validating the information by actually executing the black- 
box software. We define the term knowledge bot as a piece 
of software that has been trained using supervised machine 
learning techniques that interrogates simulation software and 
captures their underlying properties. The goal is to provide 
an accurate model of black-box physics-based simulations via 
training such a bot using artificial neural networks. 

II. Relevant Work 

The computational model of artificial neural networks has 
been widely popular in the early 1940s up until the late 
1950s for mathematics and algorithms [1], [2]. However, a new 
spark of interest emerged in the 1990s and 2000s due to its 
roots in biological relevance, where human brain computation 
and conventional digital computation differ. Artificial neural 
systems are often analogized as “physical cellular systems 
which acquire, store, and utilize experiential knowledge [3].” 
These cellular systems have been modeled by a network of 
interconnected nodes, generalizable to an input layer, one 
or more hidden layers, and an output layer. The network 
formally defined as a Multilayer Perceptron (MLP) neural 
network, where each node takes in a set of input parameters 
and computes an output value feeding to nodes further down 
the layer hierarchy as input. Eventually, at the output layer 
a value is returned as a prediction [4]. This model of the 
network acts as a foundation for numerous machine learning 
techniques using artificial neural networks. It was shown in 
Ref. [5] that multilayer feed- forward networks are a class of 
universal approximators, making the use of artificial neural 
networks sufficiently good for a large number of varying 
applications. For example, neural networks have been used 
in areas of robotics and control dynamics for the controls of 
robot manipulators and non-linear systems [6], [7] and motion 
planning [8]; in areas of computer vision for face detection 
[9], [10]; in areas of chemistry to predict secondary structures 
of proteins [11]; in areas of clinical medicine [12]; and even 
in areas of finance to predict financial markets and bankruptcy 
[13], [14], [15]. Therefore, there has been a large amount of 
research using these networks in the form of machine learning 
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Fig. 1. The role of the knowledge bot in system analysis and design process 


Fig. 2. Abstract model of bot use-case, user provided with abstraction of bot 


to capture specific patterns that can be exploited in making 
predictions, similar to humans, although in a much larger 
scale. GPU-based implementation of using artificial neural 
networks was shown in [16] to be largely successful in pattern 
recognition, allowing for massively parallel implementation. 
Various means of training neural networks to visualize and 
understand complex datasets has been demonstrated in [17]. 
However, for the purpose of our approach, we are more 
concerned with feasibility, precision, and accuracy rather than 
visualization because dimensionality reductions may remove 
essential components of our problems. 

III. Learning Physics-Based Simulations 

The knowledge bot described in Sec. I plays an important 
role in the design and system analysis process. Firstly, let’s 
assume that the knowledge bot does not exist in the process and 
we look at the main cycle of the design process illustrated in 
Fig. 1. The design process without the knowledge bot is simply 
the traditional design process of input-output relationship with 
the simulation software. 

A. Design Process without Knowledge Bot 

1) User wishes to evaluate the outcomes if several parame- 
ters are modified, so user requests the simulation of the 
system with a set of input parameters 

2) The black-box software executes the simulation and 
provides a set of output 

3) Execution is complete and output is returned to the user 

*) Repeat. Eventually the user will need to re-evaluate the 

system with another set of modified input parameters, 
thus returning to Step 1. 

Now observing this cycle, we see that the most time consuming 
event is Step 2. The physical simulation itself, may take a 
long time to execute. Furthermore, the user, if lacking expert 
intuition, may not be able to proceed with the design process 
until the outcome of the simulation is known. Additionally, 
if given sufficiently bad input parameters, the simulation may 
result in failure. However, only after the simulation is finished 
executing does the user finally discover such essential piece of 
information. The purpose of the knowledge bot is to alleviate 
such ineffectiveness by providing an instantaneous prediction 
of whether or not the simulation will execute without failure, 


such that unnecessary failing cases do not have to be simulated 
in the design process. 

B. Design Process with Knowledge Bot 

Al) User wishes to evaluate the outcomes if several param- 
eters are modified, so user requests the simulation of 
the system with a set of input parameters. However, the 
input is given to the knowledge bot. 

A2) The knowledge bot make a prediction whether simula- 
tion will be successful or failure 
A3) If knowledge bot predicts successful simulation, user 
will perform the simulation. Otherwise the input param- 
eters are adjusted and return to Step A2 
A4) Predictive output is returned to the user 
*) Repeat. Eventually the user will need to re-evaluate the 
system with another set of modified input parameters, 
thus return to Step Al. 

The role of the knowledge bot is to provide a speedy prediction 
of the outcome of simulation without executing the simulation 
itself, using the given parameters. Thus, this response is 
significantly faster than having to execute the simulation. The 
user has a prediction of the outcome of the system, thus 
the design process can proceed knowing the general outcome 
provided from the knowledge bot. The training process for a 
knowledge bot is defined as: 

LI) The knowledge bot uses a Monte Carlo sampling tech- 
nique to select candidate input parameters. The tech- 
nique uses pseudo-random numbers to select samples 
from a probability distribution. The Monte Carlo tech- 
nique has some similarity to the Latin Hypercube tech- 
nique, where samples are selected in a stratified domain 
with input probability distributions. 

L2) The software finishes executing the simulation with the 
specified parameters and provides output. 

L3) The knowledge bot uses the output to refine the learned 
model. 

*) Repeat or stop. The knowledge bot either continues to 
learn from the simulation software or user decides that 
the learned model is sufficiently accurate and stops the 
learning process. 

The purpose of the knowledge bot is to learn the underlying 
properties of physics-based simulations that govern the phys- 
ical world to help accelerate the analysis and design process. 
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Fig. 3. Underlying mechanics of the knowledge bot, abstracted away in Fig. 2 


Rather than having an analyst wait until after the simulation to 
realize that the input parameters results in a failed simulation, 
the knowledge bot is capable of providing a rapid prediction. 


IV. Knowledge Bot Mechanics 

The knowledge bot represents an abstraction of simulation 
software, where to the user, the bot is simply a trained expert 
providing accurate predictions of simulation software behav- 
iors based on user input parameters. The underlying predictive 
model, mechanics, and properties are abstracted, as illustrated 
in Fig. 2. However the knowledge bot creation encompasses 
various processes as shown in Fig. 3. The knowledge bot can 
be in a state of either training or testing. In the training state, 
the knowledge bot uses sufficiently large input data set and 
executes the simulation software. The outputs along with the 
sampled parameters are then collected in a training set. In our 
case, sufficiently large data set was denoted as 10000 instances. 
As it will be shown later, it is possible to create an accurate 
knowledge bot with less than 100 instances of simulations. 
The training set is then provided to a Java machine learning 
tool, Weka to develop the predictive model that is then fed 
back to the knowledge bot to make predictions. Similarly in 
the case of training, the bot again uses candidate samples 
input parameters, generating another sufficiently large set as 
the test set. The set is then provided to Weka, which is tested 
against the current predictive model of the knowledge bot. The 
resulting statistics are then given as output. MATLAB scripts 
is used to develop the classification results of each instance in 
the test set. We then used support vector machine algorithm 
to determine the boundary between simulation successes and 
failures. We also used MATLAB to collapse the resulting 
multidimensional classifications into two dimensions for not 
only visualization but also to describe in a single function the 
behavior of the classification boundary. 

The abstraction of the knowledge bot allows for multipur- 
pose simulation learning, as separate bots are trained to capture 
knowledge in various physics-based simulations. In fact, the 
concept of the knowledge bot can be implemented for any 
arbitrary physics-based simulations allowing for predictions of 
a wide range of physical problems and simulations. 

V. Differential Equation Classification 


Partial differential equations are generally used in physics- 
based analysis of the physical world. The proposed knowledge 
bot is capable of learning the idiosyncrasies of various nu- 
merical methods that are used to solve such partial differential 
equations as well as capturing underlying properties of com- 
plex physics-based simulations. 

Nominal classifiers were developed from a Multi-layer Per- 
ceptron neural network that is trained via the error back- 
propagation training algorithm. The problem formulation is 
based on a set of input and output vectors and a collection of 
hidden layers. The weight associated with links between layers 
are adjusted as learning rules propagate backwards from the 
output to input layer [3]. Furthermore, the learning rate for all 
sessions, unless specified otherwise is generally r = 0.3. 

Weka 1 was used for training our neural network in this 
paper; Weka is a Java-based general purpose machine learning 
software developed at the University of Waikato, New Zealand. 


Most physics-based simulations can be accurately modeled 
by a set of differential equations. Especially in the case of 
aerospace vehicle design, fluid mechanics, structural mechan- 
ics, and heat transfer are some of the most important concepts 
taken into account in vehicle systems analysis. Furthermore, 
these important concepts are governed primarily by a set of 
differential equations that share a similar form. The below 
simple partial differential equation in two dimensions is used 
to describe the classifications. 


d 2 u d 2 u d 2 u du du 

a d^ + b d^K, + c W + d^ + % 


(i) 


This can be classified into more specific families or classes of 
differential equations [18]. Given the set of coefficients a, 6, c, 
it has been shown that r = b 2 — 4 ac can be used to classify 
these classes of equations based on the sign of r as follows, 
1) parabolic differential equation, if r = 0 


1 http ://www.cs. waikato . ac . nz/ ml/weka/ 
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Classifier Performance For Each Instance 



a 


(a) Classification of each instance classified described by (1) 

* Correctly classified as hyperbolic 

■ Corrrectly classified as elliptic 

+ Correctly classified as parabolic 

* Incorrectly classified as hyperbolic 

D Incorrectly classified as elliptic 

o Incorrectly classified as parabolic 

(b) Classification legend 
Fig. 4. Classification and legend 


2) hyperbolic differential equation, if r < 0 

3) elliptic differential equation, if r > 0 

Given a simulation governed by a differential equation of the 
form described in (1), a method that classifies the equation 
into three distinctive classes and provides whether the input 
to the simulation will result in failure or success, similar to 
the acquired knowledge of a trained expert is demonstrated. 
The applications of numerical methods for these differential 
equations are treated as a piece of software with unknown 
underlying properties. The goal is to train knowledge bot to 
learn the underlying solution patterns and classify whether 
input to each of these three forms of differential equations 
will result in success or failure in numerical convergence and 
solution accuracy. 

Provided that a partial differential equation is in the form 
described in (1), it can be classified into the three classes 
described earlier. This problem of classifying the type of 
differential equation is isomorphic to classifying whether the 
solutions to a quadratic equation is undefined, real, or complex 
(however, note that the a, 6, c in the case of the quadratic 
classification corresponds to the coefficients of an equation 
in the form ax 2 + bx + c). A classifier was trained using 
10000 instances generated by Monte Carlo sampling where 
— 1 < a, 6, c < 1. The artificial neural network model used 
for training consisted of a single hidden layer with twenty-five 
sigmoid nodes. The input node of the network corresponds to 
the a, 6, c. 


When the classifier, after being trained by a set of training 
data, is given a test set to classify, it may classify each instance 
as either, 

1) r = 0, depicted as a diamond in Fig. 4(a) 

2) r > 0, depicted as a square in Fig. 4(a) 

3) r < 0, depicted as a circle in Fig. 4(a) 

The classification of the 10000 instances was 98.8% ac- 
curate using a single hidden layer network architecture with 
twenty-five sigmoid nodes. Furthermore in Fig. 4(a), an in- 
stance may either be colored solid green or outlined as red. 
The green denotes that the prediction is correctly classified. 
The red outline denotes that the prediction was incorrect. Each 
of these cases corresponds to the classification of a class of 
differential equation (parabolic, hyperbolic, or elliptic). Refer 
to Fig. 4(b) for a pictorial legend. Note that this classifier 
simply classifies a partial differential equation in the form 
described by (1) as a general class or form, whether it be 
parabolic, hyperbolic, or elliptic. It does not suggest anything 
about the equation’s stability or convergence. Simply given a 
general form of a partial differential equation, we can classify 
it as one of three classes of equations. A classifier that learns 
the convergence and stability of each of these forms will be 
discussed in following sections. 


VI. Parabolic Equation: 1D-Heat Equation 

To understand the properties of parabolic differential equa- 
tions, a distinguishable example problem was selected - the 
heat equation problem simplified to a one-dimensional space 
with initial temperature v{x) = Vo and boundary conditions 
with the temperature at zero. Thus, the physics of the problem 
assumes that the temperature change propagates from the 
boundaries, governed by the differential equation, 


dv d 2 u 

pC dt = 


( 2 ) 


It can be shown by numerical differentiation and finite dif- 
ferences (using the FTCS, forward time central-space method 
described in [22]) that the following numerical method follows, 


i >? +1 = < + 


KAt 
pcAx 2 


K+i ~ 2v" + v?_i) 


(3) 


The exact solution, however, described in [21] in the case of 
constant initial temperature v(x) = Vo, is, 


v = ^ V ^Ve-^ 1 ) 2 ^/* 2 sin { — % 1)TO (4) 


n=0 


2n 1 


In our case, we assume l = 1 and n = Kj pc. The exact 
solution is evaluated as described in (4) to the n th term where 
we decided that any terms that exceed 6 = e~ 20 is no longer 
significant for the purpose of our computations. Thus, n max 
is determined as follows, 


a/ 20/^7 (ft7T 2 £) — 1 

ft max — 


(5) 


The purpose of computing both the numerical solution and 
the exact solution (up to the n th significant term) is that we 
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Classifier Performance For Each Instance 


Q. 

5 


Ax 2 /At 


Fig. 5. Classification of each instance for parabolic lD-heat equation 


would like to treat the explicit numerical method as a black- 
box, such that it is a procedure that an individual without the 
knowledge of the heat equation or the procedure’s underlying 
numerical methods can execute and produce the temperature 
distribution after sometime t for some positions 0 < x < 1. 
Because the discretization in Eq. (3) is explicit in time, the 
algorithm will be unstable for some Ax and At combinations. 
Therefore, a classifier was trained to take the same input as the 
finite-difference procedure and to determine input validity and 
convergence of the numerical method such that an individual 
can query the classifier or software bot and determine whether 
a set of input will even execute correctly or fail on the black- 
box procedure prior to actually executing the code. The input 
parameters of the numerical method as well as the software 
bot are as follows, 

1) N x - number of discretized points in x 

Ax = 1 /(N x - 1) 

2) N t - number of discretized point in t 

At = 1 /(N t - 1) 

3) Vo - initial temperature of the substance (°C) 

4) K - conductivity of the substance 
((cal /s)/ (cm 2o C / cm)) 

5) p - density of the substance ( g/cm 3 ) 

6) c - specific heat of the substance (cal/g°C) 


* Correctly classified as success or convergent 
■ Correctly classified as failure or nonconvergent 
o Incorrectly classified as success 

□ Incorrectly classified as failure 


Fig. 6. Convergent-divergent classification legend 


The classification of 10000 numerical simulations was over 
99.2% accurate using a single layer network with twenty- 
five sigmoid nodes. It was shown by [23] that the explicit 
finite-difference numerical method has a convergence limit 
expression of a(At/Ax 2 ) < 1/2, where in this case a = 
K/ pc. Thus, the deciding term K/ pc is plotted against the 
inverse of the step ratio resulting in Ax 2 /At. The curve 
K/ pc = Ax 2 /2A t is depicted in Fig. 5. 

Note, that the term divergence is used very loosely (in this 
section and in future sections); in this context, it is more alias 
to the term non-convergence which is denoted as the result as 
one of the following, 

1) the input results in a failed execution or overflowed 
outputs 

2) the input results in divergence (or oscillation) and 
results in arbitrarily large error, e -A oc 

3) the input results in convergence but to a point that 
results in error that exceeds the maximum, 5 > £max 

VII. Hyperbolic Equation: Wave Equation 

The wave equation can be used as an identifiable example of 
the class of hyperbolic partial differential equations. The one- 
dimensional wave equation, described in [19], which governs 
the simulation of vibrations of a string or rod, can be expressed 
as the partial differential equation, 

d 2 u d 2 u 

d^ 2 = ~df 2 (6) 

where u(x,t) represents the deflection at time t of a point 
x away from the origin. The hyperbolic wave equation is 
solved numerically using FTCS; the numerical method in [19] 
is shown to be, 

< +1 = P <+ 1 + 2(1 - p)K + K -1 - < _1 ( 7 ) 

where, 



We used properties of platinum as an upper bounds for density, 
water for specific heat, and silver for conductivity for our 
training set and normalized the initial temperature, resulting 
in 0 < p < 22, 0 < iT, c, Vo < 1, and 0 < N x ,N t < 101. 
The training set for each of these networks, then is a seven- 
dimensional dataset encompassing the six input parameters and 
a convergent-divergent Boolean flag. 

After training 10000 instances, the predictions of the clas- 
sifier is illustrated in Fig. 5. The legend description is shown 
separately in Fig. 6, and teh following description will be used 
for the rest of the paper: 


The analytic solution for wave equation is, 

u(x,t) = ^[f(x + t) + f(x-t)] (9) 

where /(•) is defined to be the initial deflection, thus f(x) = 
u(x, 0). In our case of training we assumed a sinusoid where 
f(x) = sin(7rx ), since it encapsulates the boundary conditions 
u(0,t) = u(l,t) = 0. A Monte Carlo sampling of 10000 
numerical solutions was used for training based on 0 < p < 2 
and varied 10 < N x < 1001 to compute N t = 1 + 1/A t, with 
At = s/pAx . The input parameters are then simply, 
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Classifier Performance For Each Instance 


Fig. 7. Classification of each instance for hyperbolic wave equation 


Fig. 8. Classification of each instance for elliptic Laplace’s equation 


1) p - the squared ratio of step sizes in t and x 

2) N x - the number of discretized x , Ax = 1/(N X — 1) 
The trained classifier was observed and compared to the the 
actual stability criteria defined to be p < 1, which can be 
interpreted as N x < N t . 

After training 10000 instances, the predictions of the clas- 
sifier is shown in Fig. 7. The legends are identical to the 
convergence predictions of the classifiers in previous sections 
(see Fig. 6). Recall that green denotes correct classification and 
circles denote that the classifier’s prediction was convergent 
and squares correspond to divergent predictions. Note that 
the misclassifications lie along the N x = N t curve. The 
classification of the 10000 instances for the hyperbolic wave 
equation problem was over 96.2% accurate using a single layer, 
twenty-five node network architecture. 


VIII. Elliptic Equation: Laplace’s Equation 

Both Laplace’s and Poisson’s equations are important partial 
differential equations in mathematical physics and engineering 
described by [19] and are examples of elliptic equations. For 
the purpose of this paper, we used Laplace’s equation as an 
example of an elliptic partial differential equation, which takes 
the form, 

d 2 u d 2 u 

d^ + W = ( } 


The analysis of steady conduction in a rectangular region 
with prescribed boundary temperatures as described in [20] 
and it is demonstrated that a scenario with boundary conditions 
where three of the four sides as well as the interior of the region 
are initialized to zero with the last side initially T(x,y) = 
f(x), which in our case we select f{pc) = ( Au)sin(7rx ) since 
it satisfies the boundary conditions. The exact solution for 
a half-sine-wave boundary temperature distribution is indeed 
described for this particular problem, 


u(x, y ) = A u 


sinh(7ry / L) 
sinh(nH/L) 


sin(7rx/L) 


(ID 


where, A u describes the amplitude of the sinusoid and as- 
sumed to be one in this study. Furthermore, the rectangular 
region is fixed as well for simplicity, thus, H = L = 1. The 
numerical solution to this elliptic partial differential problem 
can be expressed using the successive over relaxation (SOR) 
method denoted as, 


<r=<i+ wi «i) < i2 > 

where the function L(-) denotes the residual defined as, 


l «j) = 


2 ( 1 + / 3 2 ) 


(13) 


with /3 = A | and terms Ti(-) and L 2 (-) are the central 
difference linear operators, 


L i (<j) = <+ i,j - 2u h + u 7~lj ( 14 > 

L 2 K,j) = <7+1 - 2 <i + »"j ' 1 ( 15 > 


The term uj governs how aggressively the residual plays a role 
in the relaxation of u lj- Monte Carlo techniques were used to 
sample the input parameters, which for this particular problem 
are the discretizing quantity N x ,N y and the weight for the 
residual uj. Therefore the input parameters are, 

1) N x - the number of discretized x , Ax = 1 /{N x — 1) 

2) N y - the number of discretized x , Ay = 1 /(N y — 1) 

3) uj - the residual weight 

The following input parameter domains were used in our 
demonstration, 0 < N x , N y < 101 and 0 < uj < 3. 

After training 10000 instances of the elliptic Laplace’s 
equation, the classifier was over 99.9% accurate when clas- 
sifying instances for convergence and non-convergence. The 
classification for each instance is demonstrated in Fig. 8, which 
demonstrates the predictions and whether they correspond to 
results of numerical method (see Fig. 6 for legend description). 
The network architecture of the classifier demonstrated is a 
single hidden layer, twenty-five node network. Note that the 
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Fig. 9. Classification of each instance for linearized viscous Burgers’ equation 


misclassifications and boundary between convergent and non- 
convergent lies approximately at uj = 2. 


IX. Burgers’ Equation 

The Burgers’ equation is a differential equation that de- 
scribes gas dynamics behavior. It includes a convective term 
that adds some complexity in the formulation. We chose 
linearized viscous Burgers’ equation for this study, and the 
equation is defined [23] as: 


du du d 2 u 
dt + ° dx ^ dx 2 


(16) 


The exact steady- state solution of the linearized viscous Burg- 
ers’ equation can be computed by the following, 


U = Uq 


\ £ Rl ( z,-i ) 

1 - 6~ Rl 


(17) 


where 



(IB) 


Using an explicit forward-time and central-space finite dif- 
ferences, the viscous linearized Burger’s equation can be 
discretized as: 


»r 


At 


2Ax 


U? +1 - 2u? + u?_! 

=» ( 4 a (19) 


Because the method is explicit, it may run into stability 
issues. We trained a software bot to learn about the stabil- 
ity behavior of Burgers’ equation for the following selected 
parameters: Ax, At, uq, c , and fi. 

We used a Monte Carlo technique to select the input 
parameters that are used solve the Burgers’ equation. The input 
parameters and resulting solution are used to train an artificial 
neural network to predict stability behavior of discretized 
Burgers’ equation. The parameter L = 1.0 was fixed as the 


specified domain of interest, thus the input parameters are as 
follows, 

1) N x - number of discretized x, Ax = 1/(N X — 1) 

2) N t - number of discretized t, At = 1 /(N t — 1) 

3) uo - the initial condition, initial velocity 

4) c - wave speed 

5) fi - viscosity of the fluid 

Some parameters were normalized, resulting in the following 
ranges of our input parameters 0 < uq,c,ii < 1. For 
the viscous linearized Burgers’ equation, we selected 0 < 
N x ,N t < 1001. When adjusting the network architecture the 
input nodes are labeled corresponding to the input parameters 
N x , N t , uq , c, and fi, respectively. 

The stability criteria described in [23] for the explicit finite- 
difference numerical method holds a convergence expression 
in terms of a heuristic stability region of, 

2u < Re a x < 2/u ( 20 ) 

where the mesh Reynolds number Re a x = cAx//jl and 
stability constraint v — cAt/Ax. To study the impact of Re a x 
and v on stability, we fixed N x = N t = 101 and uo = L = 1. 
We used Monte Carlo sample distributions for Re a x and v. 
The domain of our problem then lies in the following area, 
0 < Re Ax < 10 and 0 < v < 3. 

After training with 10000 instances, using a single layer with 
twenty-five sigmoid node network, the classification to predict 
stability was over 99.6% accurate. Results are shown in Fig. 9, 
which also incorporates the heuristic stability defined by Eq. 
(20) illustrated as the two blue curves. See Fig. 6 for legend 
description. Note that the heuristic stability described by Eq. 
(20) is much more conservative than capability of the numer- 
ical method as shown in Fig. 9, where the method is stable 
in regions outside the heuristic stability. The artificial neural 
network predicts the stability regions far more accurately than 
the heuristic approach. 

X. Error Behaviors in Numerical Methods 

The numerical methods presented in previous sections are 
approximate solutions to the differential equations. In order 
to evaluate and learn the accuracy behavior of these nu- 
merical methods, software hots were trained to learn and 
predict how accurate numerical solutions for a given a set 
of input parameters. The same numerical methods and exact 
solutions as presented in the previous sections were used but 
we provide the error between the numerical method and the 
exact solution as opposed to a convergence-divergence Boolean 
flag in the training data. The error is computed as the l 2 - 
norm or Euclidean norm of the difference between the two 
solutions. Furthermore, if the numerical method diverges, then 
the method may exhibit arbitrarily large error (until overflow 
occurs). Thus if divergence begins to occur, we assume that 
the error is Sm ax = 1.0, which is an upper bounds on the error 
for our training set, therefore, the absolute error is normalized 
such that 0 < e < 1. 

A numerical predictor was trained to learn the l ' 2 -norm error 
behavior for each numerical method to predict the accuracy of 
convergence for specified input parameters. In other words, 



Accuracy of Numeric Predictor Errors 




Fig. 10. Root mean squared error of the numerical predictor as network Fig. 11. Classification of each instance for the slosh propellant dynamics 

complexity increases for each numerical method convergence error predictor pendulum model problem 


the predictor acts as a convergence error predictor for each 
numerical method, predicting the error between the numerical 
solution and the exact solution. The root mean squared error 
of the predictor is illustrated in Fig. 10 as the number of 
nodes in the hidden layer increases. For this demonstration, 
a single hidden layer network with varying sigmoid nodes is 
used. It is shown that the predictors are able to predict the 
error of the numerical method within s nrn = 0.05 for methods 
in the form described in (1) of Sec. V. Burgers’ equation, 
resulted in a numerical prediction error that was slightly higher 
with e nrn = 0.07. We see that the true e nrn varies from 
method to method, and that as the number of nodes increase in 
the network architecture, the predictors plateaus at their own 
respective e nrn . The error can be further reduced by including 
additional network layers. 

XI. Slosh Propellant Dynamics 

It has been demonstrated in previous sections that artificial 
neural networks can be used to predict and classify stability 
criteria of numerical methods for computational fluid dynamics 
in addition to predicting the success and failure of simulations 
when given a set of input parameters for simulation software. 
However, the use of neural network training can be extended 
to learning other complex models. 

For large launch vehicles, a significant portion of the gross 
mass is taken up by liquid propellant, thus fuel sloshing 
dynamics plays an important role in vehicle stability [24]. 
The issue of slosh dynamics is a crucial factor for spacecraft 
and various other large flight vehicles. Sloshing dynamics may 
results in anomalies in flight with partial loss of control due to 
oscillations resulting from propellant sloshing citeslosh-motiv 
and [26]. In [28], it was shown that proportional-integral- 
derivative (PID) control systems can be designed to provide 
further stability and “avoiding adverse interaction with the 
propellant slosh and structural dynamics.” It was shown later 
that the slosh dynamics in a spacecraft can be modeled and 


simulated using a pole-cart system, in which the pole-cart 
acts as an inverted pendulum and attached to it is a second 
pendulum which acts as the slosh load [24]; the problem 
was simplified to observe the system stability in terms of the 
hinge location (the position where the second pendulum is 
attached) and the length of the second pendulum with the goal 
to establish an analytical expression that predicts the second 
pendulum’s instability. 

We used knowledge hot to classify whether parameters of 
this system may correspond to instability in the propellant 
slosh dynamics. The pendulum system is treated as a simu- 
lation with the parameters defined as, 

1) L - hinge location (position 2nd pendulum is attached) 

2) 1 2 - the length of the second pendulum 

The data given by [24] was interpolated to generate a denser 
model. Monte Carlo sampling technique was used to generate 
the training set for the artificial neural network; the sampling 
was within the domain of 0m < L < 1.0m; 0m < I 2 < 1.5m. 

After training with a set of 10000 instances, the accuracy 
of artificial neural network classification was over 99.4% 
accurate. The classification details are illustrated in Fig. 11. 
Note that the instances classified as unstable (depicted as green 
squares) lie generally in the I 2 < 0.5 region. Furthermore, 
misclassifications lie along the boundary of stable and unstable 
instances. The results demonstrate that artificial neural network 
can accurately predict complicated stability region. 

XII. Pressure Vessels Design Concepts 

The same methodologies of training can be demonstrated 
with commercial simulation software. For this demonstration, 
we observe the design of pressure vessels using NASTRAN 
[30], a commercial finite element analysis simulation pro- 
gram initially developed for NASA. The design concept of 
pressure vessels for planetary probe missions depends on 
key parameters such as the external pressure and allowable 
stress the system undergoes. In terms of systems analysis, the 
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Pressure Vessel Design Analysis 



Fig. 12. Pressure vessel design contour with a and FVdominance 

design of the pressure vessel’s mass is highly correlated to 
the maximum allowable stress a and the modulus elasticity 
E of the material. Generally, when designing the vessel, a 
material with sufficiently high a is selected in order to meet 
constraints, however, in reality, there exists instances in which 
E could dominate the design. Thus, a classifier or knowledge 
hot is used to classify and predict the instances for which the 
two governing parameters a and E dominate in contribution 
to pressure vessel design. For this particular problem we used 
NASTRAN (MSC commercial Finite-Element Analysis tool 
[29]) for simulating optimal design mass when given the 
parameters a and E. 

The mass m of the pressure vessel varies inversely with 
the root of the modulus elasticity \[E as described by [31]. 
Therefore, observing the contour plots of m\[E in the material 
domain 0 < E < 2400000, 0 < E/a < 500 demonstrated in 
Fig. 12, demonstrates that when m\fE < 450 is constant, 
corresponding to the region of E-dominance. Conversely, the 
region described by mVE > 450 then corresponds to the 
region of a-dominance. Therefore, the magnitude of , \[E was 
used as a identifier for whether the region was ^-dominant or 
cr-dominant. By cr, ^-dominances, we refer to the fact that 
either a or E is the driving parameter in terms of design and 
mass contribution. 

The training for the pressure vessel design concepts using 
NASTRAN will be discussed in a later section. But firstly, we 
will observe another type of simulation software and introduce 
a method to visually demonstrate the classifications of multi- 
dimensional results. In Sec. XIV, we will integrate this with 
the previously discussed training sets as well as the NASTRAN 
pressure vessel design classifications. 

XIII. Planetary Direct Entry Simulations 

The method of atmospheric-assisted direct entry for 
aerospace vehicles is a widely explored concept, governed by 
complex dynamics of systems of nonlinear partial differential 
equations. An analytic solution for a simplified version of the 
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Fig. 13. Classification of each instance for direct entry simulations 


direct entry problem is described in [32]. We used POST2 
(Program to Optimize Simulated Trajectories II), a point mass, 
discrete parameter targeting and optimization program. The 
goal for this simulation was to predict capture into Venus 
atmosphere, making atmospheric-assisted direct entry possible. 
Four potential input parameters that are most essential to the 
direct entry problem were observed, 

1) d - the diameter of the capsule 

2) m - the mass of the capsule 

3) v e - the entry velocity of the capsule 

4) 7 e - the entry flight path angle 

The Monte Carlo sampling techniques was used in the domains 
of, 0.5m < d < 4.0 m; 500 kg < rn < 5000 kg', 10km/ s < 
v e < 13 km/s\ and —30° < y e < 0°. 

POST2 solves a system of nonlinear partial differential 
equations of motion and generates a trajectory for the capsule 
given a set of input parameters. We determine if the input 
parameters allowed for successful direct entry by observing 
the resulting trajectory. If the trajectory resulted in a final 
position within e < 1 km of the surface of the body, then we 
consider the capsule to be captured into the atmosphere and 
as a result a successful direct entry trial. For the purpose of 
demonstration, a direct entry simulation by a 45° sphere-cone 
capsule into the atmosphere of the planet Venus was observed 
by generating 10000 instances of training data via Monte Carlo 
sampling within the specified domains for input parameters to 
interrogate the POST2 software. 

An artificial neural network was trained using the 10000 
training instances. Furthermore, we used support vector ma- 
chine (SVM) to identify the boundary between the trajectory 
cases with successful and failed direct entry parameter. The 
instances on the boundary was then used for nonlinear least 
squares regression to identify a nonlinear expression to classify 
successful direct entry cases in the following form, 

J(d,m,v e ) = am h v ( /d d 


( 21 ) 
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(a) Sec. V classifier (b) Sec. VI classifier (c) Sec. VII classifier (d) Sec. VIII classifier (e) Sec. IX classifier 

Fig. 14. Estimated boundaries of numerical method classifiers using support vector machine regression on close-perspective classification plots 


Solving for the exponents a, 6, c, and d using least squares, 
the resulting function J(-) can be used to classify the prob- 
lem of atmospheric-assisted direct entry given the parameters 

d, m,u e , and 7 e such that, 

J(d, m, v e ) = 0.0219m 0,0096 Ug ,6646 d -0,0017 (22) 

1) if 7 e < J(d, m, v e ), predicted skipped entry 

2) if 7e > J(d,m,v e ), predicted successful capture 
Furthermore, note that J(-) is defined in terms of m,v e ,d. The 
reason behind this is that y e is found to the be most influencing 
parameter in successful direct entries, thus we chose to isolate 
this in two-dimensional space resulting in J(-) defined in terms 
of the remaining parameters. 

The resulting classification on a separate set of 10000 
Monte Carlo sample input parameters is depicted in Fig. 13. 
The curve described in (22) is illustrated as the boundary 
between instances that were success direct entry parameters 
and instances that resulted in failure (skipped entries). The 
legend is similar to the previous sections with convergence 
relating to the success of the direct entry. The classification 
for the success of direct entry simulations given d, m,u e , and 
7 e using a single layer twenty-five sigmoid node network was 
over 99.6% accurate. 

XIV. Support Vector Machine Regression 

It was shown in the previous section that using the support 
vectors generated by support vector machine, combined with 
a nonlinear least squares regression, we can calculate an ana- 
lytical approximate boundary between different classifications 
(e.g., successful and unsuccessful instances). Furthermore, this 
process enabled us to reduce our multidimensional problem 
into two dimensions for more effective visualization. In this 
section we revisit the numerical method classifiers described 
in the previous sections and capture the boundary instances 
via support vector machine to compare the curve generated by 
regressions of these instances against the known analytical and 
heuristic boundaries. 

Support vector machine is a supervised machine learning 
technique that constructs one or more hyperplanes in mul- 
tidimensional space. The algorithm was originally proposed 
by Vladimir N. Vapnik [34] in 1995. The algorithm has 
been widely used in various machine learning applications, 
including text classification [35]; image retrieval [36] and 
[36]; biology and bioinformatics predicting protein subcellular 


localization [38], protein primary, secondary structures [39] 
and [40]; and bankruptcy prediction [41]. Additionally, other 
applications of the technique demonstrates its classification 
ability in terms of brain state [42] and drug use classifica- 
tion [43]. Although support vector machine can produce an 
accurate classifier for a properly- selected dataset, the tech- 
nique was chosen due to its ability to identify boundaries 
among discrete classifications. The identified boundaries can 
be used to develop an analytical expression representing the 
boundary between failure and success cases. Similar to the 
previous section, we used S VM to reduce our multidimensional 
classification results in two dimensions for more effective 
visualization. 

Support vector machine was used with previous classifiers 
to identify the boundary of success and failure (e.g., conver- 
gence or non- convergence) such that an approximate analytic 
boundary can be developed and compared to the classifier’s 
boundary. Both polynomial interpolation and nonlinear regres- 
sion were used on the classification results demonstrated in the 
previous sections. 

Let do denote the mathematically derived boundaries (e.g., 
heuristic stability boundary) and d\ denote the boundary curve 
predicted by using polynomial fit and 82 denote the boundary 
curve determined by nonlinear regression on the instances 
identified by SVM. The nonlinear regression was completed 
using the model of the form, 


Y r 


param 


= C + a(X, 


param 


r 


(23) 


where X pararn , Y pararn denote the expressions denoted in the 
x-axis and y-axis of the classification plot and the unknown 
constants ( 7 , <r, and (3 was calculated using nonlinear regres- 
sion. Furthermore, piecewise regressions and fitting was done 
on the classifier described in Fig. 14(e). Note in Fig. 14, do 
is labeled by a blue curve, d\ by a red curve and 82 by 
a black curve, as a result the following approximations for 
the numerical method classifiers in the previous sections are 
identified, 

1) Differential equation classifier, Fig. 14(a) 

do b 2 — 4ac = 0 

81 b 2 - 4ac = 0.0117 + 0.0029a 

82 b 2 - 4ac = 0.444 - 0.555a £ , 5 = 1.675F722 

2) Parabolic lD-heat equation classifier, Fig. 14(b) 

do K/ pc = 0.5(Ax 2 /A t) 
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di K/pc = 0.5131(Ax 2 /A£) - 0.0016 

62 K/pc = 0.652(Ax 2 /At) 1,237 

3) Hyperbolic wave equation classifier, Fig. 14(c) 

d 0 N t = N x 

<9i N t = 1.0019(7^) - 1.2713 

d 2 N t = 1.025 (TV*) 0 - 996 

4) Elliptic Laplace’s equation classifier, Fig. 14(d) 
do No analytic stability criteria 

<9i cj = 1.9855 - 0.0002/3 
b 2 u = 1.493 + 0.493/3 0 00036 

5) Burger’s equation classifier, Fig. 14(e) 

do Re Ax = 2 /^ 

di Re ax = 0.38z/ 4 - 6.4z/ 3 + 25. 7v 2 - 41.5z/ + 26.9 

d 2 Re Ax = 5.130z/- 1306 

do ReAx = 2^ 

di Re Ax = — 0.39z/ 2 + IMu - 0.003 

d 2 Re Ax = 1.402z/°- 855 

The boundary predictions of both polynomial and nonlinear 
regression fitting, Si,d 2 in cases 1-4, produced via regression 
on the support vector instances, were sufficiently good ap- 
proximations to the analytical boundary. The insets in Fig. 14 
demonstrate a closer view of the classification plots, demon- 
strating both the analytic boundary do as well as the predicted 
boundaries d\ and ^ 2 . 

Note that the regression for the classifier demonstrated 
in Sec. IX for Burgers’ Equation is done via a split at 
Re Ax = 1-5, resulting in two curves for the upper and 
lower constraints for the stability of the numerical method. 
The resulting equations for the upper and lower constraints 
demonstrated by polynomial interpolation (di) are given as a 
fourth order and second order equation. The approximations 
for di and &2 provide a sufficiently good border estimations, 
illustrated in Fig. 14(e) as red curves and black curves respec- 
tively. The support vector machine polynomial interpolation 
(9i, red) curves and the nonlinear regression (^ 2 , black) pro- 
vides much closer fit to the stability border than the heuristic 
stability constraints do described in (20) and illustrated in 
Fig. 14(e) as blue curves. Generally speaking via observations, 
the polynomial fitting technique producing boundary curves 
much closer to the true boundaries with well-behaved, mostly- 
linear boundaries. However, in the case of Burgers’ stability 
domain, the nonlinear regression with powers performs slightly 
more accurate as demonstrated in Fig. 14(e) due to its highly- 
nonlinear boundary. 

A. Pressure Vessel Design using NASTRAN 

Fig. 15 illustrates the classifications results using 10000 
instances for training and 10000 for testing. The resulting 
classifier, obtained by a single hidden layer network with 
twenty-five sigmoid nodes, was over 99.5% accurate. Note 
that the classification described here do not denote failure 
and success rather the regions of ^-dominance (green circles) 
and cr-dominance (green squares) are illustrated (with previous 
conventions circle corresponded to success and square to 
failure, thus one may observe the plots with ^-dominance 


UJ 
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Fig. 15. Classification of each instance for the pressure vessel simulation 


corresponding to success in mind). The red outlines represent 
incorrectly classified instances which were either identified as 
E and a but were in reality the opposite. 

The boundary equation between successful and failure cases 
was found via support vector machine and polynomial interpo- 
lation resulting in the following fourth-order equation relating 
the modulus elasticity E to the ratio E/cr at the boundary, 

Jx{E) = -0.28 E a + 1AE 3 - 2.7 E 2 + 58.9.E + 77.6 (24) 

The technique of polynomial interpolation (red) was com- 
pared to nonlinear regression (black) illustrated in Fig. 15. 
The nonlinear regression of the boundary curve in terms of 
the power equation, described in (23) for the pressure vessel 
NASTRAN simulations is, 

J 2 (E) = .031977 £ 0 - 51747 (25) 

Thus, -HE) = E /a corresponds to the instances along the 
boundary (illustrated as a red line for polynomial interpolation 
and black for nonlinear regression in Fig. 15). Furthermore, 
similar to the case described in the previous section, the 
function J(-) can be used to classify instances as either cr- 
dominant or ^-dominant generally by the following guidelines, 

1) if E/a > J(E ), predicted cr-dominant 

2) if E/a < J(E ), predicted ^-dominant 

Figure 15 shows the boundary generated by nonlinear regres- 
sion in the specified power form (black) provided a closer fit 
to the true boundary in the case of the NASTRAN simulation 
problem as opposed to the boundary generated by polynomial 
interpolation (red). 

XV. Impact of Number of Instances on Accuracy 

It should be noted that the training and testing performed in 
all previous sections with artificial neural networks used sets 
of 10000 instances. We studied the impact of training size on 
the accuracy of the learned artificial neural network model for 
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Accuracy Growth of Classifiers by Training Set Size 



Fig. 16. Accuracy growth with set size increase for mentioned classifiers 

the previous discussed problems at discretized set sizes for the 
training dataset (10, 15, 20, 50, 100, 1000, 1000). 

Figure 16 illustrates the accuracy of the classifiers as func- 
tion of training set size. With a minimum of ten training 
instances, the classifiers have a mean accuracy of 80% with 
the minimum being the hyperbolic-wave equation classifier 
at 71.77%. At 10000 training instances, the resulting neural 
network classifier achieves over 99% accuracy, with the excep- 
tion of the hyperbolic- wave equation (96.2%); the classifiers 
all exhibit similar behavior in achieving over 90% accuracy 
with less than 100 training instances and exhibits growth in 
accuracy, which corresponds to eventual convergence to all 
solutions of the black-box simulation software. Furthermore, 
some problems are presumably better suited and easier to pre- 
dict with artificial neural networks. For example, the Laplace 
equation and pressure-vessel NASTRAN classifiers were over 
90% accurate when trained with only 10 training instances. 
The legend depicted illustrates the maximum and minimum 
accuracy values between the range of 10 to 10000 instances in 
the training set. All classifiers demonstrated are achieved via 
an artificial neural network with a single hidden layer, twenty- 
five sigmoid node architecture. 

XVI. Conclusion 

We have shown that machine learning techniques can be 
used to learn the idiosyncrasies of computational simulations 
and identify regions of instability without including any ad- 
ditional information about mathematical modeling of partial 
differential equations. We have also demonstrated that the use 
of artificial neural network to learn the underlying properties 
of numerical method for solving differential equations is both 
feasible and accurate. The applications in this study included 
various types of physical simulation of fluid and heat flows 
important to systems analysis of an aerospace vehicle. 

By training classifiers using neural networks, we develop 
knowledge hots to detect simulation failures prior to execution 
and identify idiosyncrasies of unknown black-box simulation 


software. Furthermore, we demonstrated that the use of support 
vector machines helps to reduce the multidimensional and 
highly nonlinear classification datasets by obtaining simple 
expressions describing the boundary between success and 
failure instances in the dataset. 
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